W metodzie Eulera stosowanej do rozwiązywania problemów niestacjonarnych, operator różniczkowy opisujący "fizykę" symulowanego zjawiska obliczany jest w poprzedniej chwili czasowej
\( \frac{u_{t+1}-u_t}{dt} - \mathcal{L}({u_t}) = {f_{t+1}} \)
Tak więc wszystkie człony równania obliczane są na podstawie stanu w chwili poprzedniej, i wykonywana jest jedynie korekta stanu w chwili aktualnej:
\( u_{t+1} = u_t + dt * \mathcal{L}(u_t)+ dt* f_{t+1} \)
Wybieramy tak zwane funkcje testujące \( w \) i przemnażamy nasze równanie
\( \left(u_{t+1},w\right) =\left( u_t + dt * \mathcal{L}(u_t)+ dt* f_{t+1},w\right) \)
Symbol \( \left(u,w\right) \) oznacza tutaj iloczyn skalarny w przestrzeni \( L^2 \) czyli całkowanie (natomiast w następnym wzorze nawias oznacza "zwyczajny" nawias)
\( \int_{\Omega} ( u_{t+1} w ) dxdy =\int_{\Omega} ( u_t + dt * \mathcal{L}(u_t)+ dt* f_{t+1}) w dxdy \)
Zauważmy że nasz problem ma podobną strukurę jak problem projekcji bitmapy
\( \left(u_{t+1}(x,y),w(x,y)\right) =\left( BITMAPA(x,y),w(x,y)\right) \)
i tak naprawdę możemy użyć tego samego kodu, który zamiast czytać piksele bitmapy, będzie próbkował prawą stronę "piksel po pikselu".
Jak wygląda więc dokładnie nasze równanie i prawa strona?
Do aproksymacji stanu naszego systemu w danej chwili czasowej użyjemy kombinacji liniowej funkcji B-spline. W tym celu wybieramy bazę dwuwymiarowych funkcji B-spline, określając wektory węzłów w kierunku osi układu współrzędnych, na przykład dwuwymiarową bazę funkcji B-spline drugiego stopnia
\( \{B^x_{i,2}(x)B^y_{j,2}(y)\}_{i=1,...,N_x;j=1,...,N_y} u_{i,j } \)
Będą one stosowane do aproksymacji symulowanego pola skalarnego aktualnej chwili czasowej
\( u(x,y;t+1)\approx \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j} B^x_i(x)B^y_j(y) \)
i w poprzedniej chwili czasowej
\( u(x,y;t)\approx \sum_{i=1,...,N_x;j=1,...,N_y} u^t_{i,j } B^x_i(x)B^y_j(y) \)
Podobnie do testowania użyjemy funkcji bazowych B-spline:
\( w(x,y) \in \{ B^x_k(x)B^y_l(y) \}_{k=1,...,N_x;l=1,...,N_y} w^{k,l } \)
Zakładając że operator różniczkowy \( {\cal L } \) opisujący "fizykę" jest liniowy, nasze równanie wygląda teraz następująco:
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y))= \)
\( = \sum_{i=1,...,N_x;j=1,...,N_y} u^{t}_{i,j } (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y))+ \)
\( +dt * \sum_{i=1,...,N_x;j=1,...,N_y} u^{t}_{i,j } \mathcal{L}(B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y))+dt*(f_{t+1}(x,y),B^x_k(x)B^y_l(y)) \\ \forall k=1,...,N_x; l=1,...,N_y \).
Nie ustalamy tutaj konkretnego problemu, wyprowadzenia zaprezentowane tutaj dotyczą dowolnego problemu fizycznego który da się zasymulować opisaną metodą. Na przykład dla problemu transportu ciepła mamy
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y))= \)
\( = \sum_{i=1,...,N_x;j=1,...,N_y} u^{t}_{i,j } (B^x_i(x)B^y_j(y),B^x_k(x)B^y_l(y))+ \).
\( +dt *\sum_{i=1,...,N_x;j=1,...,N_y}u^{t}_{i,j }\left(\left(\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial x^2}+dt*\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial y^2 }\right),B^x_k(x)B^y_l(y)\right)+ \)
\( +(f_{t+1}(x,y),B^x_k(x)B^y_l(y)) \quad \forall k=1,...,N_x; l=1,...,N_y \)
czyli pamiętając że symbol \( \left(u,w\right) \) oznacza tutaj iloczyn skalarny w przestrzeni \( L^2 \) czyli całkowanie, mamy (w następnym wzorze nawias oznacza "zwyczajny" nawias)
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } \int_{\Omega} B^x_i(x)B^y_j(y) B^x_k(x)B^y_l(y) dxdy = \)
\( = {\color{blue} \sum_{i=1,...,N_x;j=1,...,N_y} u^{t}_{i,j } \int_{\Omega} B^x_i(x)B^y_j(y)B^x_k(x)B^y_l(y) dxdy}+ \)
\( +dt *{\color{red} \sum_{i=1,...,N_x;j=1,...,N_y}u^{t}_{i,j } \int_{\Omega } \left( \frac{\partial^2 (B^x_i(x)B^y_j(y))} {\partial x^2 }+\frac{\partial^2 (B^x_i(x)B^y_j(y)) }{\partial y^2 } \right) B^x_k(x)B^y_l(y) dxdy}+ \)
\( +dt*\int_{\Omega}(f_{t+1}(x,y)B^x_k(x)B^y_l(y))dxdy \quad \forall k=1,...,N_x; l=1,...,N_y \)
Wracając do analogii z bitmapą, mamy
\( \int_{\Omega} u_{t+1}(x,y) w(x,y) dxdy =\int_{\Omega} BITMAPA(x,y)w(x,y) dxdy \)
czyli nasza "funkcja bitmapy" dana jest
\( BITMAPA(x,y)={\color{blue} \sum_{i=1,...,N_x;j=1,...,N_y} u^{t}_{i,j } B^x_i(x)B^y_j(y) } + \)
\( +{\color{red} \sum_{i=1,...,N_x;j=1,...,N_y} u^{t}_{i,j } dt \left( \frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial x^2}+\frac{\partial^2 (B^x_i(x)B^y_j(y))}{\partial y^2 } \right) } +dt*f_{t+1}(x,y) \)
i tak naprawdę możemy użyć tego samego kodu, który zamiast czytać piksele bitmapy, będzie próbkował powyższą prawą stronę "piksel po pikselu".
Zazwczaj po prawej stronie możemy odcałkować przez części czerwony człon (zakładając że dziedzina jest na tyle duża że całka po brzegu może zostać zignorowana, ponieważ modelowane pole lub jej pochodne są tam zerowe)
\( \sum_{i=1,...,N_x;j=1,...,N_y} u^{t+1}_{i,j } \int_{\Omega} B^x_i(x)B^y_j(y) B^x_k(x)B^y_l(y) dxdy = \)
\( ={\color{blue} \sum_{i=1,...,N_x;j=1,...,N_y} u^{t}_{i,j } \int_{\Omega} B^x_i(x)B^y_j(y)B^x_k(x)B^y_l(y) dxdy}+ \)
\( -dt *{\color{red} \sum_{i=1,...,N_x;j=1,...,N_y}u^{t}_{i,j } \int_{\Omega }\left( \frac{\partial (B^x_i(x)B^y_j(y))}{\partial x}+\frac{\partial (B^x_i(x)B^y_j(y))}{\partial y } \right) \left( \frac{\partial (B^x_k(x)B^y_l(y))}{\partial x}+\frac{\partial (B^x_k(x)B^y_l(y))}{\partial y } \right) dxdy}+ \)
\( +dt*\int_{\Omega}(f_{t+1}(x,y)B^x_k(x)B^y_l(y))dxdy \quad \forall k=1,...,N_x; l=1,...,N_y \)

Zauważmy że człony (czerwony i niebieski) w których obliczamy wartości rozwiązania z chwili poprzedniej, oznaczają próbkowanie sumy wszystkich B-spline'ów przemnożonych przez współczynniki obliczone w poprzednim kroku czasowym. Całkowanie takie przeprowadza się zgodnie z rozdziałem "Sformułowanie wariacyjne a całkowanie numeryczne". W szczególności zauważmy, że całki te rozbić można na sumę całek po poszczególnych elementach siatki, i na danym elemencie siatki określonych jest jedynie \( (p+1)^2 \) funkcji bazowych B-spline, gdzie \( p \) oznacza stopień funkcji B-spine.

Symulacje metodą explicite mają ograniczenie na rozmiar kroku czasowego. Zwiększanie dokładności przestrzennej (zwiększanie rozmiaru siatki metody elementów skończonych do rozwiązania problemu projekcji prawej strony na lewą stronę) wymaga zmniejszania rozmiaru kroku czasowego, w przeciwnym wypadku symulacja zacznie zachowywać się niestabilnie (mówiąc potocznie symulacja "eksploduje"). Jest to wyrażone przez warunek Courant'a-Friedrichsa-Lewiego (CFL).
\( \frac{u dt}{dx} < C_{max } \)
gdzie \( u \) oznacza maksymalną wartość symulowanego zjawiska (tradycyjnie oznacza maksymalną prędkość dzięki temu warunek ma postać bezwymiarową co dotyczy tylko zadań z dominującą konwekcją i stałym lub zmiennym w czasie polem prędkości), \( dt \) oznacza długość kroku czasowego, \( dx \) oznacza rozmiar elementu (średnicę elementu). Warunek CFL został oryginalnie zaproponowany dla metody różnic skończonych, jednakże działa on również dla metody elementów skończonych.
Stała \( C_{max } \) nie może zostać przekroczona przez iloraz
\( \frac{u dt}{dx} \).

Jeśli tak się stanie (na przykład krok czasowy będzie za duży) wówczas symulacja zacznie zachowywać się niestabilnie (mówiąc potocznie "eksploduje") [1]

Poniżej załączam przykłady dwóch symulacji propagacji fal dla materiału spręzystego policzone metodą explicite. W pierwszym przypadku warunek CFL nie jest spełniony, w drugim przypadku jest spełniony.
Zostały one policzone kodem IGA-ADS [2]